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Abstract 

The nonlincarities found in molecular networks usually prevent 
mathematical analysis of network behaviour, which has largely been 
studied by numerical simulation. This can lead to difficult problems 
of parameter determination. However, molecular networks give rise, 
through mass-action kinetics, to polynomial dynamical systems, whose 
steady states are zeros of a set of polynomial equations. These equa- 
tions may be analysed by algebraic methods, in which parameters are 
treated as symbolic expressions whose numerical values do not have 
to be known in advance. For instance, an "invariant" of a network 
is a polynomial expression on selected state variables that vanishes 
in any steady state. Invariants have been found that encode key 
network properties and that discriminate between different network 
structures. Although invariants may be calculated by computational 
algebraic methods, such as Grobner bases, these become computa- 
tionally infeasible for biologically realistic networks. Here, we exploit 
Chemical Reaction Network Theory (CRNT) to develop an efficient 
procedure for calculating invariants that are linear combinations of 
"complexes", or the monomials coming from mass action. We show 
how this procedure can be used in proving earlier results of Horn and 
Jackson and of Shinar and Feinberg for networks of deficiency at most 



one. We then apply our method to enzyme bifunctionality, including 
the bacterial EnvZ/OmpR osmolarity regulator and the mammalian 
6-phosphofructo-2-kinase / fructose-2,6-bisphosphatase glycolytic regu- 
lator, whose networks have deficiencies up to four. We show that bi- 
functionality leads to different forms of concentration control that are 
robust to changes in initial conditions or total amounts. Finally, we 
outline a systematic procedure for using complex-linear invariants to 
analyse molecular networks of any deficiency. 



A molecular interaction network within a cell may be decomposed into ele- 
mentary biochemical reactions, such as 



where the Si are distinct chemical species. (This stoichiometry is unlikely 
but helpful for illustrative purposes.) Under mass-action kinetics, the rate 
of such a reaction is proportional to the concentrations of the substrates, 
taking stoichiometry into account. Hence, 



where x% is the concentration of species Si, Xi = [Si], and k± G M>o is 
the positive mass-action rate constant. In a biochemical network, reactions 
contribute production and consumption terms consisting of monomials like 
4Kixfx| to the rates of formation of the species in the network. This results 
in a system of ordinary differential equations (ODEs), dx/dt = f(x;n), in 
which each component rate function fi(x; k) is a polynomial in the state vari- 
ables x\, X2, ■ ■ ■ , x n € M and k%, ■ ■ ■ ,n p £ R>o are positive rate constants. 

The nonlinearities in ([2]) usually preclude mathematical analysis of the 
dynamical behaviour of such ODE systems, which are customarily studied by 
numerical simulation. This requires that the rate constants be given numer- 
ical values, which in most cases are neither known nor readily measurable. 
The resulting "parameter problem" remains a major difficulty in exploit- 
ing mathematical models, |G10| . However, the steady states of such ODEs 
are zeros of a set of polynomial equations, /i(x, k) = 0, • • • ,f n (x,n) = 0. 
Computational algebra and algebraic geometry provide powerful tools for 
studying these solutions, |CLQ97| . and these have recently been used to 
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gain new biological insights, [CDSS081 IDGHLADGlUl IMG08I IPMDSC121 
TG09a, TG09b]. The rate constants can now be treated as symbolic pa- 
rameters, whose numerical values do not need to be known in advance. The 
capability to rise above the parameter problem allows more general results 
to be obtained than can be expected from numerical simulation, [TG09bJ. 

The focus on steady states, rather than transient dynamics, is still of 
substantial interest. For instance, in time-scale separation, which has been a 
widespread method of simplification in biochemistry and molecular biology, 
a fast sub-system is assumed to be at steady state with respect to a slower 
environment and steady-state analysis is used to eliminate the internal com- 
plexity in the sub-system, |G12] , Approximate or quasi-steady states have 
also been shown to exist under various cellular conditions and can now be 
engineered in vivo, [LK99, KVBAWF04J. Finally, steady states provide the 
skeleton around which the transient dynamics unfolds, so knowledge of the 
former can be helpful for understanding the latter. 

The present paper focusses on the algebraic concept of an "invariant": a 
polynomial expression on selected state variables that is zero in any steady 
state, with the coefficients of the expression being rational expressions in 
the symbolic rate constants, [MG08j . Recall that a rational expression 
is a quotient of two polynomials; an example of such being the classical 
Michaelis-Menten constant of an enzyme, |G-B95j. (A more general defini- 
tion of an invariant allows the coefficients to include conserved quantities, 
|XG11| . but this extension is not discussed here.) Since each of the rate 
functions, fi(x;n), is zero in any steady state, the force of the definition 
comes from the restriction to "selected state variables". It is possible that, 
by performing appropriate algebraic operations on /i, • • • , f n , non-selected 
variables can be eliminated, leaving a polynomial expression on only the 
selected variables that must be zero in any steady state. 

Invariants turn out to be surprisingly useful. They have been shown 
to characterise the biochemical networks underlying multisite protein phos- 
phorylation, [MG08J, suggesting that different network architectures can be 
identified through experimental measurements at steady state. If an invari- 
ant has only a single selected variable that appears linearly, this variable 
has the same value in any steady state since it is determined solely by the 
rate constants. In particular, its value is unaffected by changes to the ini- 
tial conditions or to the total amounts of any species. This is "absolute 
concentration robustness" (ACR), as introduced in |SF1Q| . which accounts 
for experimental findings in some bacterial bifunctional enzymes, |BG03, 
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ISMMA07|fSRA09| . The mammalian Afunctional enzyme, 6-phosphofructo- 
2-kinase/fructose-2,6-bisphosphatase (PFK-2/FBPase-2), which has a more 
complex enzymatic network, also yields invariants, with implications for reg- 
ulation of glycolysis, [DCHLADG10J. The methods developed here provide 
a systematic way to analyse such bifunctional enzymes, as explained below. 

Computational algebra exploits the method of Grobner bases to provide 
an Elimination Theorem, |CLQ97j . that permits variables to be system- 
atically eliminated among the rate equations, fi,--- ,f n , [MG08J. Algo- 
rithms for calculating Grobner bases are available in general-purpose tools 
like Mathematica, Matlab and Maple and in specialised mathematical pack- 
ages such as Singular and MacaulaySj]] However, these algorithms are com- 
putationally expensive for the task at hand. They have been developed for 
general sets of polynomials and have not been optimised for those coming 
from biochemical networks. For instance, Mathematica's Grobner basis algo- 
rithm does not terminate on the network for PFK-2/FBPase-2. If invariants 
are to be exploited further, alternative approaches are needed. 

The nonlinearity in mass action comes from the pattern of substrate 
stoichiometry in ([T]) , which gives rise to the monomial nonlinearity in ([2]) . In 
the language of Chemical Reaction Network Theory (CRNT), the patterns 
of stoichiometry that appear on either side of a reaction arrow are called 
"complexes", |F79[ |G"03] , Reaction ([T]) has three species, S\,S2 and 53 and 
two complexes, 2Si + 3S2 and 45*3. If C = eiS\ + • • • + e n S n , where a £ 
Z>o are nonnegative integer stoichiometries, then the complex C gives the 
monomial, x , where x c = x^ 1 • • • x^ n . 

Aside from this nonlinearity, the defining rate equations come from linear 
processes on complexes. This observation is the starting point of CRNT 
and reveals that biochemical networks conceal much linearity behind their 
nonlinearity, [F79, G03, and see below]. This suggests the possibility of 
using fast linear methods, in preference to slow polynomial algorithms, to 
construct a subset of invariants: those that are symbolic linear combinations 
of the complex monomials, x c . As before, this definition acquires substance 
by restricting the complexes that can appear. If Ci , ■ ■ ■ ,Ck are the selected 
complexes, then a complex-linear invariant is a polynomial expression of the 
form a\x Cl + • • • + a^x k , that is zero in any steady state, where a±, • • • , Ofe 
may be rational expressions in the symbolic rate constants. 

In this paper, we examine a large class of complex-linear invariants that 



1 Available from www.singular.uni-kl.de and www.math.uiuc.edu/Macaulay2 
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we call "type 1". We determine the dimension of the space of type 1 in- 
variants (Proposition [TJ and provide a linear algorithm for calculating them 
(Theorem [TJ . We point out how invariants can be used in proving previous 
results of Horn and Jackson, |HJ72| . and the Shinar-Feinberg Theorem for 
ACR, |SF10| . We then apply the method to contrast two examples of en- 
zymatic bifunctionality, the bacterial EnvZ / OmpR osmolarity regulator and 
the mammalian PFK-2/FBPase-2 glycolytic regulator. The method is suffi- 
ciently straightforward that the invariants for networks of this kind can be 
found by manual inspection of an appropriate matrix. Finally, we outline a 
systematic procedure for analysing any network using type 1 complex-linear 
invariants. An appendix contains supplementary information. 

2 RESULTS 

2.1 Background on CRNT 

It is assumed that there are n species, S = {S\, ■ ■ ■ ,S n }, whose concentra- 
tions are xi, ■ ■ ■ , x n E M, respectively, and m complexes, C = {C±, • • • , C m }. 

c 

Complexes are regarded formally as multisets of species, C, 6 N°, where 
the value of the multiset, Cj, on the species Sj, denoted Ci(Sj) is the stoi- 
chiometry of Sj in Cj. Accordingly, Cj = Ci(Si)S\ + • • • + Ci(S n )S n . The 
reactions in the network define a directed graph on the complexes, with 
an edge Cj — > Cj whenever there is a reaction with substrate stoichiome- 
try given by Cj and product stoichiometry given by Cj. The corresponding 
mass-action rate constant gives each edge a label, Kij, treated as a positive 
symbol, «y G M>o- 

Any directed graph, G, with labels in M>o gives rise to an abstract dy- 
namics in which each edge is treated as if it were a first-order chemical 
reaction with its label as rate constant. Since the rates are all first-order, 
the dynamics are linear and may therefore be written in matrix terms as 
dy/dt = C(G).y, where y £ W 71 is a column vector, consisting of an abstract 
concentration yi at each node i of G, and C(G) isamxm matrix called the 
Laplacian matrix of G. Here, "." signifies matrix multiplication, regarding 
vectors as matrices of one row or one column. The graph Laplacian has wide 
application in biology, as described in |G12] . where additional information 
may be found. 

In CRNT, the Laplacian, C{G) : M m — > M m , provides a linear analogue 
for complexes of the nonlinear function, / : M n — > M. n , for species, in the fol- 
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lowing sense. Let : M n — > W 71 be the nonlinear function that lists the mono- 
mials for each complex, = (x , • • • , x Cm )^ . Here, ' denotes transpose. 
Let Y : W l — > W 1 be the linear function that associates to each complex, 
considered as a basis element of W 71 , its corresponding stoichiometry pat- 
tern; the i-th column of the resulting matrix is then (Cj(S'i), • • • , Cj(5 n ))^. 
With these definitions, it may be checked that f(x) = Y.C{G).'^{x) 1 for any 
x 6 M" - , as depicted in the commutative diagram in the top left of Figure [TJ 

This fundamental decomposition is due to Horn and Jackson, |HJ72] . and 
is the starting point of CRNT. They did not use the Laplacian description, 
which was introduced in [CDSS08] and exploited further in |G121 lTG09aj, 
To analyse steady states, where f(x) = 0, it is particularly useful to know 
the kernel of C(G): ker£(G) = {y G R m | C(G).y = 0}. This was first 
determined by Feinberg and Horn, |FH77[ Appendix], by a non-constructive 
method. Here, we briefly describe the constructive method introduced in 
[G12J, which shows how ker£(G) can be algorithmically calculated from G. 

This can be done in two stages, |G12| . First, if G is "strongly connected", 
so that any two distinct nodes are linked by a contiguous series of edges in the 
same direction, then dimker£(G) = 1. The Matrix- Tree Theorem provides 
an explicit construction of a basis element, pa £ M m , in terms of the spanning 
trees of G: ker£(Cr) = {pc)- The components {pc)i are polynomials in the 
symbolic labels but the details of their calculation are not needed here. If 
G is not strongly connected, it can be partitioned into its maximal strongly- 
connected sub-graphs, or "strongly connected components" (SCCs). These 
inherit from G a directed graph structure, G, in which there is an edge in G 
from SCC G u to SCC G v whenever there is an edge in G from some node in 
G u to some node in G„. G cannot have any directed cycles and so always 
has terminal SCCs, with no edges leaving them. Let these be G±, ■ ■ ■ , Gt- 
For each 1 < t < T, let p* € M. m be the vector which, for vertices of G that 
lie in Gt, agrees with the vector pc t , coming from the Matrix- Tree Theorem 
applied to Gt as an isolated graph, and, for all other vertices, j, (p t )j = 0. 
Then, the p l form a basis for ker£(G): 

ker£(G) = (/,••• , p T ) . (3) 

This gives algebraic expressions for the components of the basis vectors as 
polynomials in the symbolic labels. Note that p l may be very sparse, being 
non-zero only for vertices in the single SCC Gt- 
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2.2 Generating complex-linear invariants 

Depending on the application, invariants may be required that involve only 
certain complexes, C^,--- , Cj fc , for instance, those involving species with 
more easily measurable concentrations. Since the indices can be permuted so 
that the complexes of interest appear first in the ordering, it can be assumed 
that invariants are sought on C\,--- , Let M be the n x m matrix 
representing the linear part of the CRNT decomposition, M = Y.C(G). A 
simple way to construct a complex-linear invariant on C\ , ■ ■ ■ , is to find 
a vector, a) G R k , such that, if (a, 0)t G R m is a extended with m — k 
zeros, (a, 0) = (ai,--- ,a&, (),•■■ ,0), then (a, 0) is in the rowspan of M. 
That is, it is a linear combination of the rows of M. If x G K" is any 
steady state of the system, so that f(x) = 0, then G kerM because 

M.V(x) = Y.£(G).V(x) = f(x) = 0. Since (a, 0) is in the rowspan of M, 
(a, 0).^(x) = 0. Hence, by definition of VP, a\x Cl + • • • + akX Ck = 0, giving 
a complex-linear invariant on C±, ■ ■ ■ , 

Not all such invariants may arise in this way. For that to happen, it 
is necessary not just for (a, 0).*(x) = whenever x is a steady state but 
for (a,0).v = for all v G kerM. The relationship between kerM and 
{fy{x) | a; is a steady state} is not straightforward. To sidestep this problem, 
we focus here only on those invariants, a\x Cl + • • • +ctfcX k , in which (a, 0) is 
in the rowspan of M. We call these type 1 complex-linear invariants. Non- 
type 1 invariants do exist, as we show in the Supporting Information (SI). 
The type 1 invariants form a vector space that we abbreviate note that 
Ik depends on C±, ■ ■ ■ ,Ck and not just on k. Two basic problems are, first, 
to determine the dimension of Ik and, second, to generate its elements. 

A simple solution to the second problem is to break the matrix M into 
the n x k sub-matrix K consisting of the first k columns of M and the 
n x (m — k) sub-matrix iV consisting of the remaining m — k columns, so 
that M = K | N . Any vector 6' G W 1 which is in the left null space of N, 
b G Nl(N), so that b.N = 0, gives an (a, 0) = b.M that is in the rowspan 
of M. The assignment b — > b.M thereby defines a surjection, Ml(N) — > 
Moreover, b x .M = b 2 .M, if, and only if, (&i - b 2 ) G AA L (M) C Af L (N). 
Hence, there is an isomorphism Ik — Nl(N)/Nl(M). If X is any n x r 
matrix, dimAA^X) = n — rkX. We conclude that dim Ik = rkM — rk N, 
which yields an efficient way to determine the dimension of Ik by Gaussian 
elimination. 

In principle, this provides an automatic procedure for identifying subsets 
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of complexes with non-trivial invariants. First determine rkM and then, 
for each subset Z C {1, • • • , m}, determine the rank of the submatrix of M 
formed by those columns not in Z. If the latter is smaller than the former, 
then there are non-trivial type 1 complex-linear invariants on the complexes 
in Z. In practice, it is usually more efficient to use biological knowledge of 
the example being studied and the question being asked to narrow the choice 
of Z. We outline such a systematic procedure for finding invariants in the 
last section. 

The method above amounts to eliminating the complexes Cfc+i, • • • , C m 
by taking linear combinations of the defining rate functions, - , f n - 

This can be biologically informative because it suggests which rate func- 
tions, and, hence, which species at steady state, determine the invariant, 
IDCHLADGIO] . 

2.3 Duality and the structure of Ik 

In this section, we present an alternative procedure for calculating complex- 
linear invariants, which is based on duality and exploits the sparsity of (|3]). 
The procedure is schematically illustrated in Figure [TJ as an aid to following 
the details. 

We start, as before, with the n x m matrix, M that represents the linear 
part of the CRNT decomposition, M = Y.C{G). Let d = dimkerM and 
let B be any m x d matrix whose columns form a basis of kerM. Then, 
M.B = and the rowspan of M and the columnspan of B are dual spaces 
of each other. If G then (a, 0) is in the rowspan of M if, and only if, 
(a, 0).B = 0. If B' is the k x d sub-matrix of B consisting of the first k rows, 
then (a, 0).B = if, and only if, a.B' = 0. Hence, type 1 invariants form the 
dual space to the columns of B' . 

Proposition 1. The space 1^ of type 1 complex-linear invariants on the 
complexes C\, ■ ■ ■ , satisfies dim/^ = rkM — rkN = k — rkB' . 

Let / = rkl?'. Note that I < mm(k,d). If I = k, then dim/j. = and 
there are no type 1 invariants on C\, ■ ■ ■ , C}.. If, however, I < k, then the 
original matrix B can be simplified in two steps. First the columns. Since the 
column rank of B' is I, elementary column operations — interchange of two 
columns, multiplication of a column by a scalar, addition of one column to 
another — can be applied to the columns of B' , to bring the last d — l columns 
to zero. If exactly the same elementary column operations are applied to the 
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full matrix B, a new matrix is obtained, which we still call B, whose columns 
still form a basis for kerM. B is now in lower-triangular block form, 



B 



B' 





* 


* 



(4) 



where, as before, B' is the k x I sub-matrix consisting of the first k rows and 
I columns. 

For the rows, since the row rank of B' is still I, there are I rows of B' 
that are linearly independent. Let U C {l,--- ,k} be the corresponding 
subset of I indices and let V C {1, • • • , k} be the subset of k — I remaining 
indices. This defines a partition of the row indices of B': U n V = and 
UU V = {1, • • ■ ,k}. Let By be the I x / sub-matrix of B' consisting of the 
rows with indices in U and By be the (k — I) X I sub-matrix consisting of 
the remaining rows of B'. Using the same notation for a) £ M. k , a.B' = 
if, and only if, ajj.^B'jj) + ay.{B' v ) = 0. Since, by construction, B' v has full 
rank and is hence invertible, this may be rewritten as 

a u = -a v .(By).(B' u )- 1 . (5) 

This gives a non-redundant procedure for generating all elements of 1^ by 
choosing a v G arbitrarily and a\j £ M. 1 to satisfy (|5). The resulting 

a) £ M fc satisfy a.B 1 = and give exactly the type 1 complex-linear invariants 
on Ci, • • • ,C k . 

Using the same notation for ty{x) £ M m , the invariants themselves are 
given by ajj.^(x)u + ay.^{x)y = 0, for any steady state x £ W 1 . Substitut- 
ing ([H]) and rearranging gives av-{^{x)v — {By) .(By) -1 (x)u) = 0. Since 
ay can be chosen arbitrarily in the dual space, we conclude that 

9(x) v = {B' v ).{B' u )- l .^{x) u , (6) 

which we summarise as follows. 

Theorem 1. Each of the k — I rows of the matrix equation in ([6]) gives an 
independent type 1 complex-linear invariant on Ci, • ■ ■ , 

This procedure relies on the choice of basis elements for ker M that make 
up the columns of B and on the choice of the subset, U, of linearly indepen- 
dent rows of B'. These choices are not critical; different ones yield different 
bases for 1^. Up to linear combinations, the same invariants are found irre- 
spective of the choices. All the calculations required are linear and can be 
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readily undertaken in any computer algebra system with the rate constants 
treated as symbols. (Mathematica was used for the calculations in the SI.) 
The coefficients are then rational expressions in the symbolic rate constants. 

2.4 Haldane relationships and the Shinar-Feinberg Theorem 

Since M = Y.C(G), kerM contains the subspace ker£(G). The structure of 
the latter is known from ([3]). This should assist in the calculation of invari- 
ants, especially when ker£(G) is close to kerM. We discuss two instances of 
this, which illustrate how complex-linear invariants are related to previous 
studies. 

Define the "dynamic deficiency" of a biochemical network, 5d G N>o, to 
be the difference in dimension between the two subspaces: 5d = dimker M — 
dimker£(G), or, equivalently, 5d = dim(ker Y fl Image C{G)). This is dif- 
ferent from the "deficiency" as usually defined in CRNT, |F79[ IG03| . which 
we call the "structural deficiency", 5s G N>o- While 5q may depend on the 
values of rate constants, 5s is independent of them. However, the former is 
more convenient for our purposes. 

It is known that 5d < 5s- Furthermore, if there is only a single termi- 
nal SCC in each connected component of G, which holds for the graph in 
Figure |6]C but not for that in Figure [|/\, then 5 D = 5 S , |F79llG03] . Recall 
that a graph is connected if any two distinct nodes are linked by a path of 
contiguous edges, ignoring directions. A connected component of G is then 
a maximal connected sub-graph. Distinct connected components are totally 
disconnected, with no edges between them. 

Suppose first that 5d = and that there is a positive steady state x E 
(M>o) ?1 . Since C(G).^(x) = 0, x is a "complex-balanced" steady state, in the 
terminology of Horn and Jackson, |HJ72] , According to ([3]), the vectors p l 
provide a basis for kerM = ker£(G) and, furthermore, (p t )j ^ if, and only 
if, Cj G Gt- Choose any terminal SCC of G, which we may suppose to be G±, 
and suppose that C\ , ■ ■ ■ , are the complexes in G\ . Choose the matrix 
B so that p l is its first column and the other p l for t > 1 are assigned to 
columns arbitrarily. By construction, B is already in lower-triangular block 
form and / = 1. Setting U = {1} and V = {2, • • • , k} the k — 1 type 1 
invariants coming from (|6]) are x Ci = ((p l )i/(p 1 )i)x Cl for 2 < i < k. It is 
not difficult to see from the structure of B that these are the only type 1 
invariants. 

These invariants may be rewritten x Ci /x Cl = (p 1 )i/{p 1 )i to resemble 
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m x d k x I 



B = 


B' 


kxd elementary column 


B' 







operations 












. * . 




* 


* 



I X 



rank B' = I k x I 



B' = 




subset of I linearly- 
independent rows 







R' nonsingular 
LJ U 



(k-l) x I 



Figure 1: Schematic illustration of the dual procedure for calculating type 
1 complex-linear invariants using The fundamental decomposition of 

CRNT is shown in the commutative diagram in the top left corner, for which 
f(x) = Y.C(G).^>(x). The matrix M corresponds to the linear part of this 
decomposition, M = Y.C{G). The matrix B has columns consisting of a 
basis for kerM. The matrix B' consists of the first k rows of B. Assuming 
that rk B' = I, elementary column operations can be applied to B to bring it 
into lower-triangular block form. By a mild abuse of notation, the upper-left 
block continues to be called B' . Finally, a subset, U, of linearly independent 
rows (magenta) of B' yields the nonsingular matrix By, while the subset V 
of remaining rows (green) yields the matrix B' v , from which type 1 complex- 
linear invariants can be calculated using rol). 
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the Haldane relationships that hold between substrates and products of a 
reaction at equilibrium, |C-B95| . It follows from the construction of fr by 
the Matrix- Tree Theorem that the right-hand side of this relationship is 
determined by the rate constants, as expected for a Haldane relationship, 
|C-B95| . Horn and Jackson introduced the concept of a complex-balanced 
steady state, in part, to recover such generalised Haldane relationships for 
networks of reactions that might be in steady state but not at thermodynamic 
equilibrium, |HJ72llGT2] . 

Now suppose that 5d = 1. Then, kerM = (XiP l i''~ iP T )i where x £ 
W 71 is any vector in ker M that is not in ker C(G). Choose B to have columns 
in the same order. Suppose that there are k complexes that are not in any 
terminal SCC and that indices are chosen so that these are C\, ■ ■ ■ , Cfc. 
Then, (p t )i = for 1 < i < k and 1 < t < T, so that B is already in lower- 
triangular block form with 1 = 1. If x G (M>o) n is a positive steady state, 
then ^(x) G kerM and ^(x)i ^ for 1 < i < m. If follows that \i / 
for 1 < i < k. We may therefore choose U = {1} and V = {2, • • • , k} and 
deduce from that x Ci = (Xi/Xi) xCl f° r 2 < i < k. 

These type 1 complex-linear invariants lead to the Theorem of Shinar 
and Feinberg on ACR, |SF10| , Suppose that the structural deficiency of a 
network satisfies <Jg = 1. Suppose further that C\ and C2 are two complexes 
that are not in any terminal SCC, whose stoichiometry differs only in species 
S q . Since 5d < ^5 it must be that either 8d = or Sz> = 1. Suppose the 
former. The p l then form a basis for kerM. Because C\ is not in any non- 
terminal SCC, v\ = for any v E kerM. However, ^(x) G kerM and, since 
x G (M>o) ?1 , ^{x)\ 7^ 0. This contradiction shows that 5d = 1- It then 
follows from the invariant above that ^x q ) C2 ^ Sq ^ Cl ^ Sq,> = Xilxi- Hence, the 
steady-state concentration of S q depends only on the rate constants and not 
on the initial conditions or the total amounts and thereby exhibits ACR, 

|SFIQ] . 

2.5 Bifunctional enzymes 

The previous calculations only exploited Theorem [T] when 1 = 1. We now 
consider examples with / > 1. Details of the calculations are given in the SI. 
The examples concern enzyme bifunctionality. Enzymes are known for being 
highly specific but some exhibit multiple activities. One form of this arises 
when a protein catalyses both a forward phosphorylation — covalent addition 
of phosphate, with ATP as the donor — and its reverse dephosphorylation — 
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hydrolysis of the phosphate group. What advantage does such bifunctionality 
bring over having two separate enzymes? 

We discuss one bacterial and one mammalian example. In Escherichia 
coli, osmolarity regulation is implemented in part by the EnvZ/OmpR two- 
component system (Figure [6]A); for references, see |SMMA07] , Here, the 
sensor kinase, EnvZ, autophosphorylates on a histidine residue and catalyses 
the transfer of the phosphate group to the aspartate residue of the response 
regulator, OmpR, which then acts as an effector. Bifunctionality arises be- 
cause EnvZ, when ATP is bound, also catalyses hydrolysis of phosphorylated 
OmpR-P. 

It was suggested early on that the unusual design of the EnvZ / OmpR sys- 
tem might keep the absolute concentration of OmpR-P stable, [RS93] . This 
was later supported by experimental and theoretical analysis, |BG03j, and 
the theoretical analysis was extended to other bifunctional two-component 
systems, [SMMA07J. These ad-hoc calculations were clarified when a core 
network for EnvZ / OmpR was found to have 6s = 1 and the Shinar-Feinberg 
Theorem could be applied to confirm ACR for OmpR-P, [SF10J. Attempts 
were made to broaden the analysis by extending the core network to in- 
clude additional reactions thought to be present. For instance, EnvZ bound 
to ADP may also dephosphorylate OmpR-P. Adding these reactions to the 
core gives a network (Figure[6j3) with 5s = 2, so that Shinar-Feinberg can no 
longer be applied. However, it was shown by direct calculation in [SMMA07, 
Supplementary Information] that this network also satisfies ACR for OmpR- 
P. 

Here, we use complex-linear invariants to confirm ACR and to find a 
formula for the absolute concentration value of OmpR-P in terms of the rate 
constants. The labelled, directed graph on the complexes has thirteen nodes 
and fifteen edges (Figure |6p). Each connected component has only a single 
terminal SCC and 6d = $s = 2. We can apply Theorem [T] to systematically 
find two new invariants. 

Corollary 1. If the complexes in the reaction network in Figure are 
ordered as shown in Figure \@jp, then the space of type 1 complex-linear in- 
variants on the complexes Ci, C3, Cs, C11 has dimension 2 and the following 
are independent invariants, 

(hM^ x c 1 _ {h + h)x c 3=0 
k5X {k u + k 12 ) x {k u + k 15 ) x " U - 
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Using the expressions for the complexes in Figure [6p, it can be seen that 
x c s _ x c :i [QmpR-P] , x Cl1 = x Cl [OmpR-P] . 

Provided that [EnvZ-ATP] = x c ' A 7^ 0, the invariants can be combined and 
simplified to yield the following expression 

[OmpR-P] = 

kik 3 k 5 (k u + k 12 )(k u + fci 5 ) (7) 

hkskiok^iku + fci 5 ) + k 2 ki 3 ki 5 (k4 + k 5 )(ku + k 12 ) ' 

This confirms that, as long as there is a positive steady state, the steady-state 
concentration of OmpR-P is not affected by changes in either the amount 
of OmpR or of EnvZ. The network exhibits ACR for OmpR-P, with the 
absolute value being given in terms of the rate constants by ([7]). 

We now turn to our second example. 6-Phosphofructo-l-kinase (PFK- 
1) is one of the key regulatory enzymes in glycolysis, converting the small 
molecule fructose-6-phosphate to fructose- 1,6-bisphosphate (Figure [3^.); for 
references, see jDCHLADGlO]. In mammalian cells, the bifunctional PFK- 
2/FBPase-2 has two domains. PFK-2 has the same substrate as PFK-1 but 
produces fructose-2,6-bisphosphate. This is a terminal metabolite that is 
not consumed by other metabolic processes. Instead, it acts as an allosteric 
effector, activating PFK-1 and inhibiting fructose-l,6-bisphosphatase, the 
reverse enzyme present in gluconeogenic cells, such as hepatocytes. The 
other domain, FBPase-2, catalyses the dephosphorylation of F2,6BP and 
produces F6P. 

Biochemical studies lead to the reaction network in Figure|3p. The kinase 
domain has an ordered, sequential mechanism and the kinase and phospha- 
tase domains operate simultaneously; for more details, see |DC HLADG10| . 
The corresponding labelled, directed graph on the complexes has fourteen 
nodes and nineteen edges (Figure [IjA.). One of the connected components has 
two terminal SCCs, 5o = 4 and 5s = 5. 

Corollary 2. If the complexes in the reaction network in Figure [3]B are 
ordered as shown in Figure then the space of type 1 complex-linear in- 
variants on the complexes Ci, C2, C4, Ce, Cg, C\\ has dimension 2 and the 
following are the independent invariants, 

kix Cl - k 2 x° 2 + (k w - h)x Ce > - (kg + k u )x C8 - k 19 x Cl1 = 

k 5 x° 4 - k 8 x C(i - k u x Cs + (his ~ k 19 )x Cl1 = . 
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The second invariant in Corollary [2] was originally discovered by ad-hoc 
algebraic calculation. It is used in [DCHLADG10J to show that, if the kinase 
dominates the phosphatase, in the sense that k\s > k±g, then the steady 
state concentration of F6P is held below a level that depends only on the 
rate constants and not on the amounts of the enzymes or the substrate. 
Conversely, if the phosphatase dominates the kinase, so that kis < k±g, 
then the steady state concentration of F2,6BP is similarly constrained below 
a level that depends only on the rate constants and not on the amounts. 
Interestingly, regulation of PFK-2/FBPase-2 by phosphorylation, under the 
influence of the insulin and glucagon, causes the kinase and phosphatase 
activities to be shifted between the regimes kis > ^19 an d fcis < The 
implications of this for control of glycolysis are discussed in [DCHLADG10J. 



2.6 A systematic procedure 

The two examples discussed above had already been analysed by other meth- 
ods, so we had an idea of which invariants to expect and which subset of 
complexes to consider. For a new network, such information may not be 
available, so how can non-trivial type 1 complex-linear invariants (simply, 



"invariants") be found? The automatic procedure outlined in i2.2 can be 
used in principle but this becomes computationally infeasible when there are 
many complexes. We have found the following systematic procedure to be 
helpful on several examples. 

First determine the matrix M and from it the dual matrix B using §2.3| 
and Figure [HJ The biological context and the question being asked typically 
suggest one or more species of interest. For the initial subset of complexes, 
Z, choose all those complexes in which the species of interest have positive 
stoichiometry. Check if there are any invariants on Z using Proposition [T] If 
not, then consider any additional complexes that have at least one species in 
common with the complexes in Z. Add each of these complexes to Z in turn, 
starting with those that introduce the fewest new species and allowing the 
number of new species to increase as slowly as possible. With each addition, 
test for invariants as before. If this fails, consider adding the new complexes 
in groups, trying, as before, to minimise the number of new species that are 
introduced. 

To demonstrate this procedure, we use a modification of the EnvZ/OmpR 
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B 

EnvZ-ADP < t EnvZ 4 ► EnvZ-ATP ► EnvZ-P 
EnvZ-P + OmpR 4 ► EnvZ-P-OmpR 4 ► EnvZ + OmpR-P 
EnvZ-ATP + OmpR-P < ► EnvZ-ATP-OmpR-P ► EnvZ-ATP + OmpR 
EnvZ-ADP + OmpR-P « ► EnvZ-ADP-OmpR-P ► EnvZ-ADP + OmpR 

C, EnvZ-ADP 
C 2 EnvZ 
C 3 EnvZ-ATP 
C 4 EnvZ-P 
C 5 EnvZ-P + OmpR 
C 6 EnvZ-P-OmpR 
C 7 EnvZ + OmpR-P 
C 8 EnvZ-ATP + OmpR-P 
C 9 EnvZ-ATP-OmpR-P 
C 10 EnvZ-ATP + OmpR 

EnvZ-ADP + OmpR-P 
C 12 EnvZ-ADP-OmpR-P 
C 13 EnvZ-ADP + OmpR 

Figure 2: Two component signalling and the E. coli osmolality network. A 
Schematic of two-component phospho-transfer between a histidine residue on 
the autophosphorylating sensor kinase (light blue) and an aspartate on the 
response regulator (red) . B Extended reaction network for the EnvZ / OmpR 
two-component osmoregulator in E. coli, following [SMMA07, Supplemen- 
tary Information]. Hyphens, as in EnvZ-ATP, indicate the formation of a 
biochemical complex between the components. C Corresponding labelled, 
directed graph on the complexes, with the terminal strongly connected com- 
ponents outlined in yellow. Each connected component has only a single 
terminal SCC. D Numbering scheme for the complexes. 
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E-ATP 

E-ATP-F6P ► E + F2,6BP 
E-F2.6BP ► E + F6P 
E-ATP-F2.6BP 

E-ATP-F2,6BP ► E-ATP + F6P 
E-ATP-F6P-F2.6BP 

y E-F2.6BP + F2,6BP 

E-ATP-F6P-F2.6BP 

* E-ATP- F6P + F6P 

Figure 3: The bifunctional enzyme 6-phosphofructo-2-kinase/fructose- 
2,6-bisphosphatase (PFK-2/FBPase-2). A Schematic of the glycol- 
ysis/gluconeogenic pathway (broad gray arrows) at the step involv- 
ing 6-phosphofructo-l-kinase (PFK-1, in red), that converts fructose-6- 
phosphate (F6P) into fructose-l,6-bisphosphate (F1,6BP), and fructose-1,6- 
bisphosphatase (FBPase-1, in dark blue), that catalyses the opposing reac- 
tion in gluconeogenic tissues. PFK-2/FBPase-2 (light blue) operates bifunc- 
tionally to produce and consume fructose-2,6-bisphosphate (F2,6BP), which 
allosterically regulates PFK-1 and FBPase-1, as shown. B The correspond- 
ing reaction network. 



B E «=► 

E-ATP + F6P 4 ► 

E + F2,6BP < ► 

E-ATP + F2.6BP < ► 
E-F2.6BP < ► 
E-ATP-F6P + F2,6BP < » 



E-ATP-F2.6BP + F6P < ► 
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B C 1 


E 


c 2 


E-ATP 


C 3 


E-ATP + F6P 


C 4 


E-ATP-F6P 


C 5 


E + F2,6BP 


C 6 


E-F2,6BP 


C 7 


E + F6P 


C 8 


E-ATP-F2.6BP 


C 9 


E-ATP + F2.6BP 


C 10 


E-ATP-F6P + F2,6BP 


C 11 


E-ATP-F6P-F2.6BP 


C12 


E-ATP-F2.6BP + F6P 


C 13 


E-F2,6BP + F2.6BP 


C 14 


E-ATP-F6P + F6P 



Figure 4: PFK-2/FBPase-2, as in Figure [3j in terms of complexes. A La- 
belled, directed graph on the complexes, with the terminal SCCs outlined in 
yellow. The last connected component has two terminal SCCs. B Numbering 
scheme for the complexes. 
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network in Figure [3^3 . We add a single new reaction 



OmpR-P — ► OmpR 

for spontaneous (non-catalysed) dephosphorylation of OmpR-P. The phos- 
pho-aspartate bond in response regulators is labile and may be spontaneously 
hydrolysed, so this new reaction is biochemically plausible. The labelled, 
directed graph for the modified network in Figure [5]A. has fifteen nodes and 
sixteen edges. Each connected component still has only a single terminal 
SCC, like the graph in Figure [6p, but now <5d = 5s = 3. The matrices M 
and B are provided in the SI, along with other details of the calculation. 

The biological context suggests that the active state of the response reg- 
ulator, OmpR-P, is of most interest. Following the procedure above and 
using the table in Figure leads to Z = {CV, C\\, C14} as an initial 
subset of complexes. It can be readily checked by inspection of B that the 
corresponding rows yield a submatrix of full rank 4, so that Proposition [T] 
tell us that there are no non-trivial invariants on Z. Among the remaining 
complexes, C\, C2 and C3 each involve only species that are already present 
among the complexes in Z. Adding each to Z in turn, it can be checked 
that each of the subsets Z U {Ci}, Z U {C2} and Z L) {C3} have a space of 
invariants of dimension 1, which respectively yield the following non-trivial 
invariants, 

k ^ x ° 14 - (45fe) xCl + x ° H + (ASfe) xCl1 = 

hex^ - x Ci + X C S + (k^\ X C U = o (8) 



k l6 x c ^ - k 5 x c * + (^_) xC* + (^) x c ^ = . 

These all have a similar form, due to the common subset Z. The absence of 
C-j in these invariants could have been inferred directly from the pattern of 
entries in B (SI). 

A non-trivial invariant does not necessarily provide helpful biological 
insights. This depends crucially on the context and the question being stud- 
ied. For instance, assuming that we have a positive steady-state, the second 
invariant in ([8]) may be rewritten in terms of steady-state species concentra- 
tions as 



[EnvZ] 



[OmpR-P] = -. , V " 4+fc5/ - r . (9) 

ho + (^fc) [EnvZ-ATP] + (^gu) [EnvZ-ADP] 
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Equation ([9j establishes a steady-state relationship between the activated 
response regulator, OmpR-P, and the sensor kinase, EnvZ, under nucleotide 
loading. This may be useful depending on the available experimental data. It 
does suggest that OmpR-P no longer exhibits ACR, as it did for the network 
in Figure [6^5, since its steady-state level appears to depend on the amount 
of EnvZ present. However, not much more can be said just from ([9]). 

The situation is different for the first and third invariants in They 
lead to a similar equation for [OmpR-P] as in (^J but with the numerators 
on the right hand side given by, respectively, 

/, ' | /, ' :; ''''' » [EnvZ-ADP] and k 5 [EnvZ- ATP] . 



k 2 {k4 + k 5 



Because the same species now appears in both the numerator and the de- 
nominator, a simple comparison and cancellation, yields the inequalities 



[OmpR-P] < < 



fci k 3 k 5 (fci4 + k 15 ) 

k 2 (fc* + h) h 3 ki 5 

(fen + fe 12 ) 
kb ~i — I • 



The strictness of the inequality comes from the assumption that the steady 
state is positive. We see that the activated response regulator has two upper 
bounds that are robust: they depend only on the rate constants and not 
on the initial conditions or the total amounts of either EnvZ or OmpR. An 



interesting aspect of (10) is the absence of parameters k§, ■ ■ ■ ,kg, that relate 
to the phosphorylation of OmpR. The other reactions contribute factors to 
the bounds that can be biochemically interpreted. Recall that the catalytic 
efficiency of an enzyme is the ratio of its catalytic rate to its Michaelis- 
Menten constant, k cat /KM, |C-B95j . The factor ^13^15/(^14 + fcis) in the 
first bound is the catalytic efficiency of OmpR-P dephosphorylation by EnvZ- 
ADP, while the factor ^10^12/(^11 + ^12) m the second bound is the catalytic 
efficiency of OmpR-P dephosphorylation by EnvZ- ATP. The balance between 
these redundant dephosphorylation routes will influence which of the two 
bounds is the tighter and this balance is further modulated by the efficiency 
of nucleotide binding to EnvZ. The first bound is modulated by the factor 
^3^5/(^4 + &s)> which may be treated as an effective catalytic efficiency for 
EnvZ phosphorylation by ATP (it has units of sec -1 rather than M _1 sec~ 1 ), 
and the factor k\/k 2 , which is the equilibrium constant for ADP binding. 
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The second bound is only modulated by k§, which is the catalytic rate for 
EnvZ phosphorylation. 

The existence of these robust bounds also strongly suggests that OmpR- 
P does not satisfy ACR. Indeed, it can be shown by algebraic calculation 
that [OmpR-P] can take different values in different steady states (SI). In 
contrast to the second invariant in (J8|, the first and third invariants yield 
interesting and unexpected biological insights. Further non-trivial invariants 
can be sought using the procedure above and we leave this for the reader. 

We note two further points of interest. First, the addition of a single 
new reaction to a network can markedly change its behaviour from exhibit- 
ing ACR to only having robust bounds. This is a feature of biochemical 
networks under mass- action kinetics. It raises difficult problems of interpre- 
tation because there is always the possibility that the actual cellular network 
may include reactions that have been missed in a model. Very little work 
has been done on this difficult question. Invariants may provide a way to 
study it: perhaps certain invariants can be shown to remain "invariant" when 
a network is enlarged in a particular way. Second, a theme is emerging from 
the examples considered here. Despite the differences in network structures 
between Figures [HJ [3] and |5j the bifunctionality in each case serves to limit 
the steady state concentration of a substrate form, either absolutely, as in 
Figure [HJ or relative to some robust upper bound, as in Figures [3] and [5] We 
speculate that this may be a design principle of those bifunctional enzymes 
that catalyse forward and reverse modifications. There are other forms of 
bifunctionality, such as enzymes that catalyse successive steps in a metabolic 
pathway, and preliminary studies suggest that these behave very differently. 
If modification bifunctionality did evolve to implement concentration control, 
markedly different network structures seem to have converged upon it. 

3 DISCUSSION 

The nonlinearity of molecular networks makes it impossible to solve their 
dynamical behaviour in closed form. Their analysis has therefore relied on 
numerical integration and simulation, for which the biochemical details and 
the numerical values of all parameters must be specified in advance. This 
has made it difficult, if not impossible, to "see the wood for the trees" and 
to discern general principles within the overwhelming molecular complexity 
of cellular processes. Invariants are part of a new approach in which net- 
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EnvZ-ADP 


c„ 
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EnvZ 


c, 
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EnvZ-ATP 


C 4 


EnvZ-P 


C 5 


EnvZ-P + OmpR 


c 6 


EnvZ-P-OmpR 


C 7 


EnvZ + OmpR-P 


C 8 


EnvZ-ATP + OmpR-P 


C 9 


EnvZ-ATP-OmpR-P 


C 10 


EnvZ-ATP + OmpR 


C 1, 


EnvZ-ADP + OmpR-P 


C 12 


EnvZ-ADP-OmpR-P 


c 13 


EnvZ-ADP + OmpR 


C 14 


OmpR-P 


C 15 


OmpR 



► 15 



Figure 5: Complexes for the modified EnvZ/OmpR network discussed in 



£2.6, in which a single reaction is added to the network in Figure [6J3. A 
Labelled, directed graph on the complexes, with the terminal SCCs outlined 
in yellow. Only the connected component at the bottom differs from the 
graph in Figure [6p. B Numbering scheme, with two additional complexes, 
C14 and C15, beyond those in Figure [6t 



work behaviour at steady state can be analysed with the parameters treated 
symbolically. There are now several examples, drawn from different biolog- 
ical contexts, in which the invariants summarise the essential steady-state 
properties of the network. The key to wide exploitation of this method is 
that it should be readily applicable to realistic networks. In principle, Grob- 
ner bases allow any invariant to be calculated but this is computationally 
infeasible in practice. Complex-linear invariants form a limited subset of all 
invariants but, as shown here, they have biological significance and can be 
efficiently calculated for realistic networks. Our work clarifies previous re- 
sults and provides a new tool for symbolic, steady-state analysis of molecular 
networks. 
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5 Appendix 



The symbolic linear algebra calculations that follow may be readily under- 
taken in any computer algebra system. We used Mathematica, for which we 
wrote custom functions that compute the relevant matrices automatically 
from a description of the reaction network. 



5.1 Invariants that are not of type 1 

Consider the hypothetical reaction network in Figure [6j\. While such chem- 
istry is unlikely, it illustrates the mathematical issues. The network has three 
species, S%,S2, S3 and nine complexes C\, ■ ■ ■ , C9, ordered as in Figure [Hp. 
The ODEs are 

dxi 

—7- = k\xi 
at 

^ = k 4 xix 3 - k 2 x\ (11) 

= k 5 xix 2 - k 3 xl . 
With the given ordering, the matrix M = Y.C{G) is 


fc 4 
-k 3 k 5 

Focussing on the complexes C±, C3, C5, C7, Cg, Proposition [T] shows that the 
space of type 1 complex-linear invariants has dimension three. However, it 




is easy to see from (11) that the only steady state of the network is when 
xi = X2 = x.3 = 0. Hence, for any values of o, b, c,d,e S K, the polynomial 
expression 

ax Cl + bx Cs + cx°' 5 + dx° 7 + ex° 8 

always vanishes in any steady state. Hence, the space of complex-linear 
invariants on C\, C3, C5, C7, Cs has dimension five. 
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Figure 6: Network with complex- linear invariants that are not of type 1. 
A Hypothetical reaction network. B Labelled, directed graph on the com- 
plexes, with the terminal strongly-connected components outlined in yellow. 
C Numbering scheme for the complexes. 
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5.2 Corollary [T] 

The nine species and thirteen complexes in the EnvZ/OmpR network in 
Figure 1 are ordered as follows. 



Si 


EnvZ-P-OmpR 


C x 


5*8 


EnvZ-ADP 


s 2 


EnvZ-ATP-OmpR-P 


c 2 


1S4 


EnvZ 


S3 


EnvZ-ADP-OmpR-P 


c 3 


s 7 


EnvZ-ATP 


s^ 


EnvZ 


c 4 


s 9 


EnvZ-P 


s 5 


OmpR 


c 5 


5*9 + S*5 


EnvZ-P + OmpR 


s e 


OmpR-P 
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EnvZ-ATP 
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EnvZ + OmpR-P 


s 8 


EnvZ-ADP 


c 8 


5*7 + ^6 
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EnvZ-P 
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EnvZ-ATP + OmpR 
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EnvZ-ADP + OmpR-P 






C12 


S3 


EnvZ-ADP-OmpR-P 






C13 


5*8 + S5 


EnvZ-ADP + OmpR 



With this ordering and with the rate constants as in Figure 1C, the matrix 
M = Y.C(G) is 
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J 








fc 5 





— fc 6 


k? 
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A basis for the kernel of M can then be calculated to make up the columns 
of a matrix B. 



h 













V 1 



fcifc 3 fc 5 
(k 4 +k 5 )k 15 
k r j k^ 

fcl5 
k 5 







k 2 (k 4 +k 5 )k 12 
kik 3 k 5 
(k 4 +k 5 )k 12 

fcig 

fc5 


n 
u 


n 
u 


n 
u 


(k 7 +k$)k 15 
kekg 

fcis 

kg 






(fc 7 +fc 8 )fci2 
fcl2 

















kio 








1 





1 





^14+^15 
fel3 








1 























k&k% 
kg 

ka 
1 















o\ 




1 










V 



Corollary [T] fo cusses on the complexes C\, C3, C$, C\\, so that k = 4. These 
are not the first four complexes in the ordering, as was assumed for conve- 
nience in §2,5| We can imagine that the columns of M and the rows of B 
have been permuted so that these complexes are now the first in the order- 
ing but we will not bother to write out these new matrices. We note that 
columns 2 and 4 of B have non-zero entries in the relevant four rows, while 
the remaining columns have zero entries. We can undertake elementary col- 
umn operations on B, as described in the paper (in fact, only interchange 
of columns is required), to bring B into lower-triangular block form. The 
resulting 4x2 sub-matrix, B' 7 in Equation Q, is then given by 

/ k 2 (k 4 +k 5 )k 15 



V 



kik 3 k 5 

kis 
k 5 



fcl4+fcl5 
kl3 



k 2 (k 4 +k 5 )k 12 \ 
kik 3 k 5 > 

ki2 
k B 
fcll+fcl2 







/ 



The columns of this are linearly independent, so that rkB' = 2. It follows 
from Proposition [T] that the dimension of the space of type 1 complex- linear 
invariants on Ci, C3, Cg, Cn is 2, as claimed. To generate the invariants, we 
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note that rows 2 and 3 of B' are linearly independent, so that we can take 
U = {2, 3} and V = {1, 4}. Then 



B' 



V 




B' 



fcl3 



k2(k4+k 5 )ki2 v 
fclfc3fc5 





Since ^f(x)i 



x 



C*3 rj.Cs 



s )t and ^{x)v = {x 1 ,x ix y , the two linearly inde- 



pendent type 1 complex-linear invariants may be read off from Equation ^ , 



v Ci 



' k 2 (k 4 +k 5 ) 
kik 3 

fc5(fcl4+fcl5) 
fcl3fci 5 



fcipfcl2(fcl4+fcl5) 
kl3k 15 (kn+k 12 ) 




to yield the expressions in Paper Corollary [TJ as claimed. 
5.3 Corollary g 

The eight species and fourteen complexes of the PFK-2/FBPase-2 network 
in Paper Figures 2 and 3 are ordered as follows. 



Si 


F2,6BP 


Ci 


s 5 


E 


s 2 


F6P 


c 2 


s 7 


E-ATP 


s 3 


E-ATP-F6P 


c 3 


5*7 + 5*2 


E-ATP + F6P 


s 4 


E-ATP-F6P-F2,6BP 




Ss 


E-ATP-F6P 


s 5 


E 


c 5 


S 5 + S! 


E + F2,6BP 


s e 


E-F2,6BP 


c 6 


Se 


E-F2,6BP 


s 7 


E-ATP 


c 7 


s 5 + s 2 


E + F6P 


s 8 


E-ATP-F2,6BP 


c 8 


Ss 


E-ATP-F2,6BP 






c 9 


S 7 + Si 


E-ATP + F2,6BP 






Cio 




E-ATP-F6P + F2,6BP 






C u 


S4. 


E-ATP-F6P-F2,6BP 






C\2 


Ss + S2 


E-ATP-F2,6BP + F6P 






C*13 


5*6 + Si 


E-F2,6BP + F2,6BP 






C14 


s 3 + s 2 


E-ATP-F6P + F6P 



With this ordering and with the rate constants in Paper Figure 3A, the 
matrix M = Y.C(G) is 









fc 5 


— fc 6 


k 7 





fcl3 


-fcl2 


-fcl4 


fcl5 + fcl8 














-fc 3 


fc4 





*>8 





fcll 








fcl7 + fcl9 


-fcl6 











fc 3 


— fc4 — /C5 

















-fcl4 


fcl5 + fe 19 



































fcl4 


— fcl5 — fcl7 — fclS — fcl9 


fcl6 








-fcl 


k 2 


fc 5 


— fcg 


k 7 + k s 



































kQ 


— k 7 — fc 8 — fc 10 













fcl8 











fcl 


— k 2 —k 3 


k-1 











fcll + fcl3 


-fcl2 





























fcio 





— fcg — fcn — fci3 


fcl2 





fcl7 


-fclS 
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and a matrix B, whose columns form a basis for the kernel of M, can then 
be calculated as 



fc 8~ fc 10 




(-k R +k ln )k 12 


fc 8 fc 9 + fc 10 fc ll 





hi k 10 


H 


k l k 10 


k 1 k 10 



















-k 5 k 10 + (k 4 +k 5 )k 8 


- fc 4 fc 18 + ( fc 4 + fc 5) fc 19 


(k 4 -\-k 5 )k 8 k 12 


(fc 4 + fc 5 )(fc s fc9+fc 10 fc 11 ) 





k 3 k 5 k 10 


fc 3 fc 5 


fc 3 fc 5 fe 10 


fe 3 fe 5 fc 10 




*=8 


-k ls + k lg 


fc 8 fc 12 


fc 8 fc 9 + fc 10 fc ll 





fc 5 fc 1Q 


k 5 


fc 5 fc in 


fc 5 fc 10 




k r + k s + k 10 


fc 18 


(fc 7 + fc 8 + fc 10 )/s 12 


(fc 7 + fc 8 )A; 9 





k 6 k 10 


fc 6 


k e k io 


k 6 k 10 




1 





_ k 12 


fc 9 





k lQ 




k 10 


















1 











1 











1 


fc ll + fc 13 











k 12 




1_ 


fe 15 + fc 18 + fc 19 











fcl4 


fc 14 









1 











1 


k L7 











^16 


k l6 







































1 













/ 



Corollary [2] focusses on the complexes C\, C2, C4, Cg, Cs, Cn, so that k = 
6. As before, we can imagine that the columns of M and the rows of B 
have been permuted to make these complexes first in the ordering. Only 
columns 3, 4, 5, 6,8 of B have non-zero entries in the relevant rows and, 
when restricted to these rows, column 5 is a scalar multiple of column 3. As 
before, we can interchange columns to bring B into lower-triangular block 
form, with the resulting 6x4 sub-matrix, B', in Equation Q given by 





fcifcio 


kw 
fcl 


fcgfcg+fciofcll 

fcifcio 


ki 1 













1 




^8 


-fcig+fcl9 


fcgfcg+fciofcll 







k s k 10 


k 5 


fcsfclO 




1 

kw 





_feg_ 

fcio 













1 





V 





1 





0/ 



with B' evidently of full rank 4. It follows from Paper Proposition 1 that 
the space of type 1 complex-linear invariants on C±, C2, C4, Cq, Cs, Cn has 
dimension two, as claimed. To generate the invariants, we can take U = 
{2, 4, 5, 6} and V = {1, 3} so that 



B f , 



u 



(0 








1\ 


1 

fcio 





kio 











1 





\o 


1 





0/ 



' k$-k 10 
fcifcio 
fcs 

v ^5^10 



fcl9 
fcl 

-fcl 8 +fcl9 



fcsfcg+fciofcn fca 
fcifcio fcl 

fcsfc9+fclQfcll 

fcsfcio 
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B'jj is evidently non-singular. The two linear-independent type 1 complex- 
linear invariants can then be read off from Equation Q, 



fx 
X'i 



&2 (fcs — fcio) (fcg+fcll) fcig 
fe'l fe'l &1 fcl 

kg kii (— ^18+^19) 

^5 ^5 ^5 







\x Cll J 



to yield the expressions in Paper Corollary |2j as claimed. 
5.4 Example of the systematic procedure 

The example in § |2.6| is a modification of that in Corollary [T] As shown 
in Figure 5 it has nine species and fifteen complexes which are ordered as 
follows. 



Si 


EnvZ-P-OmpR 


Ci 


5*8 


EnvZ-ADP 


s 2 


EnvZ-ATP-OmpR-P 


c 2 


1S4 


EnvZ 


s 3 


EnvZ-ADP-OmpR-P 




S 7 


EnvZ-ATP 


s 4 


EnvZ 


c 4 


5*9 


EnvZ-P 


s 5 


OmpR 


c 5 


5*9 + 5*5 


EnvZ-P + OmpR 


s 6 


OmpR-P 


c 6 


Si 


EnvZ-P-OmpR 


s 7 


EnvZ-ATP 


cv 


S4 + S*6 


EnvZ + OmpR-P 


Ss 


EnvZ-ADP 




S7 + Sq 


EnvZ-ATP + OmpR-P 


s 9 


EnvZ-P 


c 9 


S2 


EnvZ-ATP- OmpR-P 








s 7 + s 5 


EnvZ-ATP + OmpR 






C u 


Ss + Sq 


EnvZ-ADP + OmpR-P 






C\2 


S3 


EnvZ-ADP-OmpR-P 






Cl3 


Sg + S5 


EnvZ-ADP + OmpR 






Cl4 


s 6 


OmpR-P 






Cl5 


s 5 


OmpR 



With this ordering and with the rate constants shown in Figure 5A, the 
matrix M extends that for Corollary [T] with a 9 x 2 block on the right: 



/ 














kti 


-kr - k a 


kn 












































feio 


-fen - fci2 

















































fcl3 


-&14 - k 15 










ki 


— k2 — kz 


ki 








ktt 


— kg 







































fey 








A' 12 








kiB 





fel6 



















fe8 


— kg 


-feio 


fcn 





-ki3 


ki4 





-&16 







k:i 


-ki - kz, 














— kio 


fcll + fcl2 

















V 


-hi 


k 2 


























-ki3 


kl4 + fcl5 














ks 





— ks 


k 7 
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A matrix B, whose columns form a basis for kerM, has the following form. 



/n 

/ o 


U 


k2 (k4,+ks)kia 
k-ik-sks 
(^4+^5)^16 

fel6 

k B 






n 
U 


^2(^4+^5)^15 
fcifc3&5 

k 3 k 5 
fcis 
fes 






n 
U 


^2(^4+^5)^12 
kiksks 

(^'4 + ^5)^12 
^3^5 
ki2 

k s 






n 
U 


\ 



n 
U 























1 






(k7+k$)ki6 
k(>k$ 

fcl6 
kg 






(fc7+fc8)fcl5 
k$k$ 
fcis 

fc8 






(k\7+ks)ki2 
kekg 
ki2 
kg 


kjktj 
k&k$ 

kji 

ks 
























1 




















fcii+fcia 
feio 























1 




















1 




















fcl4+fcl5 
fcl3 























1 




















1 




















1 




















V 




















0/ 



Consider the initial set of complexes, Z = {C7, C§, Cn, C14}, suggested by 
the procedure in the Paper. The corresponding rows of S have only a single 
non-zero entry, each of which occurs in a distinct column, so that the four 
rows constitute a sub-matrix of full rank. It follows from Proposition [T] 
that the space of type 1 complex-linear invariants is empty. Next consider 
adding in turn each of the complexes C\, C2 and C3 to Z. Consider first 
ZL){Ci}. The rows of B may be permuted so that the rows corresponding to 
complexes C\, Cg,Cn, C14, C-j are placed first and in that order. Permuting 
the columns of B so that the columns 2, 4, 6, 7 are placed first yields a matrix 
in lower-triangular block form, with the top left block, B', given by 





k2(k4+k 5 )k 16 
kik 3 k 5 





k 2 (k4+k 5 )k 15 
k\k 3 k 5 



^14+^15 
fcl3 


k 2 {k4+ks)ki2 
kik^ks 

fell+fcl2 

feio 












1 













V 











1 


/ 
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It is evident that B' is in block diagonal form, with the last row and column, 
corresponding to complex Cy, forming an identify matrix. It is then easy 
to see from Equation ^ that Cy cannot appear in any invariant. In other 
words, it is sufficient to work with Cs,Cu,Ci4, when adding C\. Before 
doing this, note that the C\, C2 and C3 rows of B may be written in the 
form 

fcis ki2 0) 



where 



aj(0 feie 

k 2 {k 4 + k 5 ) 
kifak 5 



a 2 



(fc 4 + k 5 ) 
k 3 k 5 



"3 



1 

k 5 



(12) 



Accordingly, we can do all three calculations at once by omitting Cy and 
re-writing B' in the form 



B' 



( aik 16 



V 1 



aik 15 


fcl4+fcl5 
fel3 





Q-ikyi \ 



fcll+fcl2 

fcio 






/ 



with Oj given by (12) for i = 1,2,3. It is evident that ikB' = 3 and we can 



choose U = {1, 2, 3} and V = {4}, so that 



B\ 



V 



It follows that 



/ aikiQ 


V 



{B'u 





fcl3 



V 



1 






atiki2 

fcu+fcig 
fcio 





fcipA;i2 

(fcn+fci2)fci6 



fcio 

fell+fel2 



B' V = {1 0). 



(fei4+fci 5 )fci6 
fcl3 
^14+^15 





Note that on only appears in a single entry. Using Paper Equation (6), with 
^{x)jj = (x Ci , x Cs , x Cl1 )^ for i = 1,2, 3, and ^f(x)v = (x Cl4 ), we recover the 
invariants in Equation (18]). 



5.5 Failure of ACR for the example in § 2.6 



The steady states of the example just discussed can be algebraically sim- 
plified as follows. The differential equations governing the system can be 
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obtained from the matrix M in §5. 4 1 above by using the fundamental decom- 



position of CRNT, dx/dt = M.^(x). With the notation in the Table in j pT4 
this gives 



dSi 

~dt 

dS 2 

dt 
dSs 

dt 
dSi 
~dt 
dS 5 

dt 
dS 6 
~dt 
dS 7 
~dt 
dS s 

dt 
dS 9 

~dt 



k 6 S 5 S 9 - (k 7 + k 8 )Si + kgSiSe 
k w S 6 S 7 - (fcn + fci 2 )5 2 

kl 3 S e S$ - (fcl4 + h 5 )S 3 

kiS s 



(13) 
(14) 
(15) 
(16) 
(17) 



(k 2 + k 3 )Si + kiSr + k s S± - k 9 SiS 6 
k 7 S\ - k 6 S 5 S 9 + k 12 S 2 + k 15 S 3 + k 16 S e 
k 8 Si - kgS^Se - k w S 6 S 7 + hiS 2 ~ ki 3 S e S$ + k u S 3 - ki 6 S 6 (18) 
k 3 Si - (fc 4 + h)S 7 - k 10 S 6 S 7 + (feu + k 12 )S 2 (19) 
k 2 Si - kiS$ - k 13 S 6 S 8 + (ku + h 5 )S 3 (20) 



k 5 S 7 - k 6 S 5 S 9 + k 7 Si 



(21) 



It can be checked that these equations (or, equivalently, the matrix M from 
which they are derived) satisfy two conservation laws, corresponding to the 
total amounts of sensor and response regulator, 



Si + 5*2 + S3 + S4 + S 7 + S$ + Sq = K\ 
S 1 + S 2 + S 3 + S 5 + S 6 = K 2 , 



(22) 



where K\, K 2 G I^>o are constants determined by the initial conditions. 

The steady states are obtained by setting the right hand sides of equations 
(13)-(21) to zero. The nine variables may be partitioned into two subsets, 
{Si, S2, S3, 5*7, Ss} and {S4, 55, Sq, Sg}, in such a way that the variables in 
the first subset can be written in terms of the variables in the second subset. 



Using equation ( 13 ) , 



Si 



k 6 



k 7 + k 8 



S5S9 + 



k 7 + 



S^Sq . 



Combining equations (14) and (19), we get 



s 7 



ki + fc 5 



S, 



(23) 



(24) 
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and, similarly, combining equations (15) and (20) we get 

S 8 ~- 



fa 



fa 



Sa . 



Using equation ( 14 ) together with ( 24 ) we get 

ho \ ( fa 



So 



hi + k 



12 



fa + fa 



S^Sq 



and, similarly, using equations (15) and (25) we get 
S 3 = 



ku + ki5 



fa 



fa 



S^Sq . 



(25) 



(26) 



(27) 



Equations (23)-(27) describe {S\, S2, S3, Sj, Ss} in terms of {S4, S5, 5*6, Sg}. 



If we now substitute in equation (17) for Si, S2 and S3 using equations 



(23), (26) and (27), respectively, and simplify, we get 



(aS 4 + fa 6 )S 6 



fa + fa 



S5S9 



(28) 



where a depends only on the rate constants and can be conveniently repre- 
sented as a = (3 + 7, where 



fak 



7K9 



fa + fa 



1 



fank 



10«12 



fai + fa 2 



fa 



fa + fa 



ki3fa 5 \ (fa 
kiA + fa 5 J \fa 



If we also substitute in equation (21) for Si and S7 using equations (23) and 



( 24 ) , respectively, and simplify, we get 
fafa 



fa + fa 



S4 + PSaSq 



fak, 



6«8 



fa + fa 



S5S9 



(29) 



Combining equations ( 28 ) and ( 29 ) we can express Sg as a rational function 
of £4, 

s& = (fa + fa)( 7 s 4 + k 16 ) ■ (30) 



Finally, substituting this expression into equation ( 28 ) , we can express S5 as 
a rational function of S4 and Sg, 



Sb 



fafajfa + fa)(aS4 + fcigjiSj 
fafa(fa + fa)(^S4 + fa e )Sg 



(31) 



35 



If S4 and Sg are given arbitrary positive values, the values of all the other 



variables are determined by equations (31), (30) and (23) to (27). It can be 



checked that these values form a positive steady-state of the system. The 
free quantities 64 and Sg, in terms of which the other variables have been 
parameterised, implicitly determine the values of the conserved quantities K\ 
and K2 in equation (22). It can be seen from equation (30) that Sq, which 
is the concentration of the activated response regulator, 56 = [OmpR-P], 
varies with the choice of S4 and, hence, does not exhibit ACR. 



It also follows from equation (30) that 



[OmpR-P] = S 6 < 



(fc 4 + k 5 )j 



(32) 



which shows, as deduced in £2.6, that OmpR-P has a robust upper bound. 



The bound in ( 32 ) is half the harmonic mean of the robust bounds in Equa- 



tion 10 and is therefore tighter than either of the latter. This is evident from 



equation ( 30 ) , where it is clear that Sq asymptotically approaches the bound 
in (32) as S4 increases. Hence, (32) is the best possible bound on Sq. 
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